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Abstract 



We p resent an extension of the multi-moment advection scheme ( iMinoshima et al. 



20111 . J. Comput. Phys.) to the three-dimensional case, for full electromag- 
netic Vlasov simulations of magnetized plasma. The scheme treats not only 
point values of a profile but also its zeroth to second order piecewise moments 
as dependent variables, and advances them on the basis of their governing 
equations. Similar to the two-dimensional scheme, the three-dimensional 
scheme can accurately solve the solid body rotation problem of a gaussian 
profile with little numerical dispersion or diffusion. This is a very impor- 
tant property for Vlasov simulations of magnetized plasma. We apply the 
scheme to electromagnetic Vlasov simulations. Propagation of linear waves 
and nonlinear evolution of the electron temperature anisotropy instability 
are successfully simulated with a good accuracy of the energy conservation. 

Keywords: Advection equation, Conservative form, Multi- moment, Vlasov 
simulations, Magnetized plasma 



1. Introduction 

The kinematics of collisionless plasma has been studied in a wide va- 
riety of fields, such as in laboratory plasma physics, space physics, and 
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astrophysics. Evolution of collisionless plasma and self-consistent electro- 
magnetic fields is fully described by the Vlasov-Maxwell (or Vlasov-Poisson) 
equations. Thanks to recent development in computational technology, self- 
consistent numerical simulations of collisionless plasma have been successfully 
performed from the first-principle Vlasov-Maxwell system of equations. 

One of the numerical simulation methods for collisionless plasma is so- 
called Vlasov simulation, in which the Vlasov equation is directly discretized 
on grid points in ph ase space. Compared t o the most popular Particle-In- 
Cell (PIC) method (IBirdsall and Langdonl . Il99ll ). the Vlasov simulation is 
free from the statistical noise inherent to the PIC method. This advantage 
can allow us to study in detail such as wave-particle interaction, particle 
acceleration, and thermal transport processes, in which a high energy tail 
in the velocity distribution function plays an important role. On the other 
hand, the Vlasov simulation requires a highly accurate scheme for the advec- 
tion equation in multidimensions, to preserve characteristics of the Vlasov 
equation (i.e., the Liouville theorem) as much as possible. It also requires 
larger computational cost than the PIC method. 

A number of advection schemes have been proposed for the application to 
the Vlasov simulat ion thus far (e.g.JCheng and Knorrill976uNakamura and Yabei . 



g 

19991 : iFilbet et all l200ll ; iMangenev et alXT2002l ; ICrouseilles et al.l . l2009h . Al 



though the schemes have been succeeded especially in applying to the elec- 
trostatic Vlasov-Poisson simulation, the application to the electromagnetic 
Vlasov simulation of magnetized plasma is still limited, mainly owing to the 
difficulty in solving the gyro motion around the magnetic field line (solid 
bo dy rotation in velocity sp ace) . 



f 

Minoshima et al.l (120 111 ) have proposed a new numerical scheme for the 



advection equation, specifically designed to solve the Vlasov equation in mag- 
netized plasma. The scheme treats not only point values of a profile but also 
its zeroth to second order piecewise moments as dependent variables, and 
advances them on the basis of their governing equations, for better conser- 
vation of the information entropy and reducing numerical diffusion. In the 
paper, we have presented one- and two-dimensional schemes, and have shown 
their quite high capabilities. Especially, the two-dimensional scheme can ac- 
curately solve the solid body rotation problem of a gaussian profile with little 
numerical dispersion or diffusion. These schemes have been successfully ap- 
plied to electrostatic and electromagnetic Vlasov simulations. 

For the application of the electromagnetic Vlasov simulation to a wide 
variety of magnetized plasma phenomena, however, it is necessary to de- 
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velop a three-dimensional scheme to treat the full three-dimensional velocity 
space. In this paper, we present an extension of our previous schemes to the 
three-dimensional case. Similar to the one- and two-dimensional schemes, 
the three-dimensional scheme treats the zeroth to second order piecewise 
moments as well as point values of a profile as dependent variables. Details 
of the scheme are described in Section [2J Benchmark tests of the scheme 
are presented in Section The application of the scheme to electromag- 
netic Vlasov simulations is presented in Section HI Finally, we summarize 
the paper in Section [5j 



2. Three-dimensional multi-moment advection scheme (MMA3D) 

We consider the time evolution of a three-dimensional profile f(x, y, z, t) 
and its zeroth to second order moments in the x, y, and z directions defined 

as 



M° = jjj fdV (= M° = M° y = M° z ) , (1) 
= ^ JJJ x m fdV, (m = 1,2) , (2) 

K = —All y m f dV > (ro = l,2), (3) 



ml 




MT = ^ / // z m fdV, (m = l,2), (4) 

where dV = dxdydz. The conservative advection equation of / and governing 
equations of the moments are written as 

af at at at /e« dv ai»\ 

-L +u— +v— +w— = - ( + + /, (5) 

dt dx dy dz \ dx dy dz J 

aj & + l dx iLUI ufdyd *) + I dy i (// vfdzdx ) + hi (// wfdxdy ) = °- (6) 

9M, 



01 



(m - 1) 
dM. 



(m - 1) 



+ / dx — — f — — ff ufdydz | -f- / dy — — ( f f vx m jdzdx \ -f- f dz f / / wx m fdxdy\ 

J dx \ m! J J J J dy \m\ J J J J dz \m\ J J J 

~TT /// uxm ~ 1 f dV ' (™ = 1 > 2 ). ( 7 ) 

+ / dx h (h. // u y mfd y dz ) + / ^ (^r // vfdzdx ) + / d ^ // ^ m / d ^) 

—— JJJ vy m ~ l fdV, (m = l,2), (8) 



if + J dx jx- (i // « m '*"«) + / // + / // 

7 fff wzm ~ lfdV - ("1=1.2). (9) 

(m — 1)! J J J 
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where u, v, and w are the velocity component in the x, y, and z directions. 
Here, the conservative advection equation of / is cast into the advective form. 
Eqs. ©-Q are obtained by multiplying Eq. (jSJ) by x m /m\, y m /m\, or z m /m\, 
and then integrating over space. Hereafter, we assume du/dx = dv/dy = 
dw/dz = 0, because we are concerned with the Vlasov equation. We use 
vector forms M m = (M™, M™, M™ ) , x = (x, y, z), and u = (u, v, w). 

To solve a set of these equations, the three-dimensional MMA scheme 
treats eight dependent variables; the point value of the profile fij t k and the 
piecewise moments, 

j r*k+i rvj+i pxi+i 
M»i/2j+iM+i/2 = — } / / (^ = 0,1,2), (10) 

m - J z h Jyj J Xi 

and constructs a piecewise interpolation for / in a cell with a quadratic 
polynominal, 

3 3 3 

Fi,j,k (x, !/. z ) = ^EE vii\C ull x t i,j,k {x - Xi)^ 1 (y - yjT" 1 (z - z k ) v ~ x , (11) 

v=\ fj,=l A=l 

which gives an interpolation function for M m as 

G%, k (x,y,z) = —J Z n X x fm F i)jjk (x',y',z')dV 

m - JziJyjJxi 
3 3 3 / A™(x,Xi) \ 

= EEE W<v>vi) 

„=i^=ia=i \ A™(z,z k ) J 

xCvn\;i,j,k ~ Xi) X (y - yjY [z - z k ) v , (m = 0, 1, 2) (12) 
where G m = (G™, G™, G™), G° x = G° y = G° z = G°, and 
(x, Xi) = 1, 

A\(x,Xi) = (Aa; + x^ j (A + 1) , ( 13 ) 
A\ (x,Xi) = |A(^i) x 2 + Ax . x + x 2| /{(A + 1)(A + 2)}> 

To determine the coefficients C UfJi \-ij±, we use the variables at the upwind 
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position as constraints, 



where 



■^i,j,k ip^ii Vji %k) 
■^i,j,k v^iupi Vji %k) 
■^i,j,k y^ii Vjupi %k) 
Fi,j,k \%ii Vji Zkup) 
Fi,j,k (%iupi Vjupi Zk) 
Fi,j,k y^ii Vjupi Zkup) 
i^iupi Vji Zkup) 
■^i,j,k {p^iupi Vjupi %kup) 
(*i,j,k (-Eiupi Vjupi z kup) 



xup 

icell 

jup 

jcell 

kup 

kcell 



fi,j,ki 

fiup,j,ki 

fi,jup,ki 

fi,j,kupi 

fiup,jup,ki 

fi,jup,kupi 

fiup,j,kupi 

fiup,jup,kupi 

sgn (Ci,j,k) sgn (rji^k) 



(14) 



sgn {6. 



icell, j cell, kcelli 



i + sgn (Q,j,k) i 
i + sgn (Q dtk ) /2, 
j + sgn (r)t,j,k) i 
j + sgn (rj iJjk ) /2, 
k + sgn (9i,j,k) , 
k + sgn (9i,j, k ) /2, 



0,1,2) 



(15) 



are the position of the upwind grid and cell in the x, y, and z directions, 
sgn (C) stands for the sign of (, and (0,j,kiVi,j,ki9i,j,k) is the distance of the 



upwind departure position relative to 
order accuracy, 



J. ; 



Vji Zk) 



determined with a second 




-Ui,j, fc At + [(« ■ V) u] 



At 2 



i,j,k 2 



Eq. (Tl4|) is obviously insufficient to determine the twenty-seven coeffi- 
cients C UfJj x-ij^- Then, we additionally introduce line-integrated variables 
I = i}xi lyi h) defined as 



2fe + l 



r x i+i rvj+i 

lx;i+l/2,j,k = / fdx, l y -i,j+l/2,k = / fdy, l Z ;i,j,k+l/2 = / fdz, (16) 

Jx t Jyj J z k 
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and use their upwind values as constraints, 






f T 

J^x;i,j,k 


i^iupi Vji %k) 


= S£ 


? n (Cm 


k) lx;icell,j,ki 




L X ;i,j,k 


{p^iupi Vjupi %k) 


= S£ 


? n (Cm 


k) lx;icell,jup,ki 






{p^iupi Vji Zkup) 


= 


? n (Cm' 


k) ' x;icell ,j ,kupi 






i%iv/pi Vjupi Zkup) 


= 


? n (0,j 


k) lx;icell,jup,kup 




Ly;i,j,k 


\p^ii Vjupi ^k) 


= si 


5 n (Vij 


k) ly;i,jcell,ki 


< 


Ly;i,j,k 


Vjupi %kup) 


= Si 


5 n (Vi,j 


k) ly;i,jcell,kupi 
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y%iupi Vjupi %k) 


= Si 


5 n (Vi,j 


k) ly;iup,jcell,ki 




Ly;i,j,k 


y%iupi Vjupi Zkup) 


= Si 


5 n (Vi,j 


k) ly;iup,jcell,kup 




Lz;i,j,k 


K %i ; Vj ; Zkup) 


= Si 


? n {hi 


k) I z;i,j ,kcelli 




Lz;i,j,k 


p^iupi Vji Zkup) 


= Si 


? n (#M 


k) I z;iup,j ,kcelh 




L Z ;i,j,k 


sp^ii Vjupi Zkup) 


= Si 


? n ihj 


k) ^z;i,jup,kcelli 




, Lz;i,j,k 


J^iupi Vjupi %kup) 


= Si 




k) lz;iup,jup,kcell 



(17) 



where 

Lx\i,j,k (x, V, z) = / Fi,j,k{x',y,z)dx' 

J Xi 

3 3 3 

= ^^^v^C vli x.i d . k {x ~ x.f {y-yjY" 1 (z- zkf 1 , (18) 

u=l n=l A=l 

and Ly^jfi, Lg-^j^ are given likewise. Consequently, the coefficients are ex- 



plicitly determined, which are listed in Appendix A 

To save memory cost of our Vlasov simulation code, we do not treat 
the line-integrated variables as dependent variables. Therefore, they should 
be constructed from known variables. We find that a construction tech- 



nique used in the Weighted ENO scheme (IJiang and Shul . 119961 1 works well 
for problems we are concerned. The line- integrated variable Z X ;i+i/2 is con- 
structed from four point values (fi-i, fi, fi+i, fi+2) as 



x;i+l/2 



Ax [a L (-/*_! + 8/, + 5/ i+ i) + a R (5/, + 8f i+1 - f i+2 )) 



12 (a L + a R ) 



(19) 



&L,R = (IS l ,r + e) 



where we use p = 2 and e = 10 6 (same as in Jiang and Shu ( 19961 )). and a 



uniform grid spacing is assumed. Simulation results are not sensitive to the 
choice of p, e. ISl,r is a smoothness measurement of the interpolation func- 
tion qi,R (x) (second order polynominal) on the left- and right-side stencils, 
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[i — 1, i, i + 1) and (i,i + l,i + 2), denned as 

'dq L ,R 



IS 



L.R 



Ax 



d 2 q L ,R 



Ax- 



<9x 5 



Qr {x 
giving 

IS L 

is R 



dx 

fi-l — 2fi + fi+1 

2Ax 2 l ~ ~" ' 2Ax 

fi — 2/j+l + fi+2 , N 2 . fi+2 — f 

— (x-x i+1 ) + 



dx. 



>2 fi+1 — fi-l / _ \ , f 
X X, J ~r n A l^X X, J ~~r 



2Ax 2 



2Ax 



(x - x i+ i) + /i+i, 



(/i-i - Vi + /i+l) (5/i-l - 16/, + ll/ m ) - fU> 



6 
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(/i - Vi+i + /< +2 ) (H/i - 16/j+i + 5f l+2 ) + {fi+2 ~ fi 



6 4 
Eq. ( IT9l) gives a fourth-order central interpolation when = or. 

2.1. Time integration for the advection 

Let us consider the time integration of the variables. For the Vlasov 
equation in magnetized plasma, we consider two problems of the advection 
with constant velocity and the solid body rotation. When the velocity is 
constant in space, we employ the semi-Lagrangian method. This method has 
been also applied to the one- and two-dimensional schemes. The advection 
phase of Eqs. ©-© is calculated as 



n+l f. . , 

n+ls|0 

ice.ll ,j cell ,k cell 



nF i,i,k 0* + C, Vj + v, Zfc + Q) , 

sgn (C) sgn {rj) sgn (9) 



(20) 



f z kup pVjup f x iup 

/ / / 

-J I I n Fiup,j,kdV 

J zu-hQJ +rw x; 



Zk+9J yj+r\J Xi up 
z kup fyjup~^~V f x i 



+ / » 

+ / / n F lJ:kup dV 

r z kup rvjup+v r^iup+C 
+ / / / 
f z ku P + fVjup+v r 

+ 

J z kup J Vjup J x 

+ / / n F, 



d\ ' 



iup,jup,k 



up,j,kup 



dv 



fZkup+8 rvjup+i rxiup+C 

+ y j 

J z kup J Vjup J X 

sgn (C) sgn (t)) sgn (0) 



"M 



icell ,j cell ,kcell 



iup,jup,kup 



dV 



(21) 



z kup rVjup r x iup~b~£ 



r z kup rvjup r^iup 

Jzk+eJyj+vJxi+C 

r-kup rvjup r*iu 

J z k +eJ vj+rfl x iup 
rzkup rvjup+v rxiup 

Jz k +eJy Jllp Jxi+C 
^ r-kup+e rvj up r^iup 

J z kup Jyj+vxi+C 
r*kup rvjup+v r*iup- 

' zk+e-'yjup J^i-up 
+ / / 

Jz kup J Vjup ■'SB<4 

r z kup+B rvjup rxiup- 
+ / / 

J z k-u.p J Vj+V J Xiup 

r z ku P + e ryjup+v r^i 

+ j y 



iup,j,k ^ 

1 jp JT/ 

r i,jup,k u,v 

l jp JT/ 

" i,j ,kup u,v 



^ " r iup,jup,k tLV 



G " ^i,jup,kupd\^ 



■* zup,j ,kup u ' v 

H 

t" 1 . ™ p. , rfV 

x lZip,JZip,K up u ' v 



(m = l,2), (22) 



where the left-superscript indicates the number of time steps and the asterisk 
means that the variables are at the intermediate step. Here, we consider the 
case of the CFL number \uAt/Ax\, \vAt/Ay\, \ wAt/Az\ < 1 for simplicity. 
Integrations on the right-hand side of Eqs. (j2ip and (I22p can be exactly 
calculated by using Eq. ( lT2"j) . Next, we advance the non-advection phase of 
Eqs. ©- " 



as 



lvl i+l/2,j+l/2,fc+l/2 

n+1 ji*-2 

JW i+l/2,j+l/2,fc+l/2 



J "i+l/2,j + l/2,*;+l/2 "+" Jw i+l/2j+l/2,fc+l/2 



uAt, 



(23) 



i+l/2j+l/2,fc+l/2 
*- M i+l/2,j + l/2,fe+l/2 + " +1 M i+l/2,j + l/2,k+l/'2~^~ J u ^ t - ( 24 ) 



2.2. Time integration for the solid body rotation 

When the velocity varies in space (including the solid body rotation prob - 
lem), we employ a time integration method proposed by lli and Xiaol (120071 ). 
in which cell-integrated values are advanced by the finite volume method 
with the Runge-Kutta time integration, whereas point values are advanced 
by the semi-Lagrangian method. This method has been also applied to the 
two-dimensional scheme. 
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In the three-dimensional solid body rotation problem, the velocity is gen- 



erally given as (u, v,w) 



ZUJy, zu x — XUJ Z , XUJ. 



where U) 



(u x , Uy, u z ) is the angular velocity. With an arbitrary u>, the rotation axis 
passes through simulation cells with an arbitrary angle, causing the velocity 
reversal within a single cell. This is an unfavorable situation for conservative- 
form upwind schemes. To avoid this situation, we split it into two phases, 
(u, v,w) = (yuj z , zu x ,xuj y ) and (—zu y ,—xuj z ,—yu x ), and then alternately 
advance. 

At the first phase, Eqs. ()6])-(j9]) are approximated into the following finite 
volume formulation, 



d M ° + l/2J+l/2,fc+l/2 

dt 



-Uj+l/2^zd x 
— z k+l/2 u; xdy 

-x i+1 / 2 oj y d z 



fdydz 



fdzdx 



fdxdy 



j+l/2J+l/2,fe+l/2 



i+l/2j+l/2,fc+l/2 



(25) 



fi i\/r m 

OIV1 i+l/2,j+l/2,k+l/2 

at 



o., 



ml 



i+l/2 j+l/2,fe+l/2 

x m fdydz 



z k +i/2^x r r r x mj dzdx 



ml 



+ 



^ 2/i+i/2W 2 M^ + 1 1/2)i+1/2jfc+1/2 

z fc+l/2 W i^+l /2,j+l/2,fc+l 12 
\ X i+l/2^yM z . i+1/2J+1/2M1/2 J 



i+l/2,j+l/2,k+l/2 

i+l/2j+l/2,fe+l/2 

i+l/2j+l/2,fc+l/2 
\ 

, (m = 1,2126) 



where d x f i+ i/2j,k = fi+i,j,k — fi,j,k- Eq. fl25l) guarantees the conservation of 
mass. Area-integrated variables of / appearing on the right-hand side of Eqs. 
( |25|) and ( 1261) are constructed from the interpolation function such as 



x m fdxdy 



sgn (Ci,j,fc) sgn (i7i,j,fc) 



icell,jcelLk 



Vjup 



Uj 



x m F itj ,k (x,y,z k )dxdy 



s g n (0,i,fc)sgn (?7ij,fc) 



p=l A=l 



, a;,-) Gi^i^fcAa^Ay", (m = 0, 1, 2) , 



where Ax = x iup - x { and Ay = y jup - y. y 
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For stable calculation, the time integration of Eqs. (I2"5"j) and (|2l)j) is im- 
plemented with the third-ord er TVD Runge-Kutta method (IShu and Osherl . 



19881 ; iGottlieb and Shul . Il998[ ). Intermediate values of / at each stage, which 



are necessary to calculate the coefficients of the interpolation function, are 
approximated by solving the equation of the characteristics with the Runge- 
Kutta method, 

d ( &,j,k \ ( {yj + Vi,j,k}u z \ ( d,j,k \ 

jj: I Vi,j,k 1 = 1 {zk + Oi,j,k} u x 1,1 \,j,k J = 0, 

fi,j,k Fi,j,k Ci,j,ki Uj ^i,i,ki @i,j,k) j (^7) 

where the left-superscript s = 0, 1, 2, 3 denotes the Runge-Kutta stage. The 
point value is advanced by = 3 fi,j,k- 

In addition, we should calculate intermediate values of I at each stage, 
which are also necessary to calculate the coefficients. Integrating Eq. (151) 
over x, the governing equation of l x is approximated as 

dl x f , d . dl x dl x 

— + / dx— (uf) + v— + w— = 0. 
at J ax ay az 

The second term is advanced by the finite volume method, and third and 
fourth terms are advanced by the semi-Lagrangian method. The solutions at 
stages are approximated as 

l X ,icell,j,k Sgn (^ijjjfc) L x -i j k (^Xi U p, t/j -\- T]icell,j,k: %k @icell,j,k) 



2 



— VjU z ddd x ° ficell,j,k, 
l X ,icell,j,k Sgn (Cij'jfe) I J x\i ) j,h i^iupi Vj ^]icell,j,ki %k @icell,j,k) 

9x ( f + 1 f)icell,j,k 

-yjUJ z At J —, 

lx,icell,j,k Sgn (Ci,j',fc) Lx>iJ t k {%iupi Vj ~\~ 1]icell,j,k^ %k @ieell,j,k) 



-yjUj z At 



«.( / + V + 4 2 /W 



6 

and l y ,l z are approximated likewise. 

Consequently, we can calculate the coefficients at each stage, and then 
advance the moments with the Runge-Kutta method, which is implemented 
as follows, 

l M m = "M m + H( n /, n i, n M m )At, 
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2 M m = ^ n M m + ~ { 1 M m + R( 1 f,H, 1 M m )At} , 

*M m = - n M m + - { 2 M m + R( 2 f, 2 l, 2 M m )At} , (m = 0,1, 2) ,(28) 
3 3 

where R stands for the right-hand side of Eqs. (1231) and ( 1231) . 

The second phase is advanced in a similar way. The time integration of 
the whole system is carried out with three steps; a half time step at the first 
phase ( n /, n M) -»■ (*/, *M), a full time step at the second phase (*/, *M) -> 
(**/, **M), and then a half time step at the first phase (**/, **M) -> ( n+1 /, n+: 



3. Benchmark tests 



As a benchmark test, we simulate the long time solid body rotation and 
advection problem, 



d f r/ \ , df n 

- + { (,-, o )x.}.^ = 0, 



of a gaussian profile, 
/ 0, y,z,t = 0) = exp 



2<r= 



2o 3 



2a? 



(29) 



(30) 



This equation describes the rotation around (x,y,z) = (xo,yo,Zo). To solve 
the equation, we split it into the rotation and advection phases, and advance 
them as follows; the advection with a half time step, the rotation with a 
full time step (this includes three steps), and then the advection with a 
half time step. The angular velocity is (u x ,u y ,co z ) = (l/y/6, l/v3, l/v2). 
The simulation domain is [—1, 1] with 32 grid points in each direction. The 
open boundary condition is employed where constant incoming fluxes are 
assumed while outgoing fluxes are perfectly lost. The time step is 27r/750. 
The simulation runs till hundred rotation period s. We compare the results 



with the CIP-CSL2 scheme ( ITakizawa et all . 120021 ). Note that both the MMA 
and CIP-CSL2 schemes treat eight dependent variables in three dimension. 
Fig. [TJ shows the results for a symmetric gaussian profile with a x = a y = 
0.2, x = y = z = (without the advection). Compared to the 



a 



CIP-CSL2 scheme (c), the MMA scheme (b) completely preserves the profile 
even after a hundred of rotations. From these simulation runs with different 
grid sizes, we evaluate the order of accuracy of the schemes. Fig. [2] shows 
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the error £\ . k \f(x, y, z, t) — f(x,y,z,0)\/N as a function of the grid size 
(N is the number of grid points). Both schemes show nearly the third order 
accuracy in space (dashed line). The error of the MMA scheme (triangles) 
is ~ 10 -L5 times smaller than the CIP-CSL2 scheme (diamonds). With the 
finest grid size, the accuracy of the MMA scheme is slightly worse than the 
third order. Currently, we cannot identify its cause. 

Fig. [3] shows the results for an asymmetric gaussian profile with <j x = 
0.15, o y = 0.2, a z = 0.25, Xo = yo = z = (without the advection). By 
fitting the profile with the gaussian function, Fig. H] shows the temporal 
variation of the standard deviation (a x , a y , a z ). While the CIP-CSL2 scheme 
shows the rapid increase due to numerical diffusion, the MMA scheme keeps 
the standard deviation with small errors. The smallest deviation a x (a) 
slightly increases, whereas others (b,c) decrease. However, their average (d) 
is kept constant. 

Fig. [5] shows the results for a symmetric gaussian profile with a x = a y = 
a z = 0.2, xq = 0.2, y = 0.15, z = —0.1. The MMA scheme provides a better 
solution with keeping (a x , <r y , a z ) and (x , y , z ) constant, indicating that the 
scheme can accurately solve the electric field (E x B) drift motion with little 
numerical dispersion or heating. Hence, the three-dimensional MMA scheme 
is a very suitable method for long time Vlasov simulations of magnetized 
plasma. 



4. Electromagnetic Vlasov simulations 

We apply the three-dimensional MMA scheme to electromagnetic Vlasov- 
Maxwell simulations. The one-dimensional electromagnetic Vlasov-Maxwell 
system of equations is written as 

d f- +Vi ^ + ^( E + ^). d fl = , (s = f) , e) , pi) 



dt dx m s \ c J dv 

E OB r 

— = cVxB- 47Tj, — = -cVxE, 3=J2^J v f* dv > ( 32 ) 

s=p,e 

where E(x) and B(x) are the electric and magnetic fields, j(x) is the current 
density, c is the speed of light, q s is the charge, m s is the mass, f s (v,x) is 
the phase space distribution function, and the subscript s denotes particle 
species (p for protons and e for electrons). Although configuration space is 



12 



assumed one dimension, full three-dimensional velocity space and electro- 
magnetic fields are treated. 

In the simulation, we treat sixteen dependent variables for both electrons 
and protons; point values of the distribution function, piecewise moments in 
the velocity space, and their cell-integrated values in the configuration space, 

fi,j,k,l f (Vx;ii Vy',ji ^ z\ki %l) i 

I f^z-k+l (""yj + l rv x ;i+l 

MT + l/2, 3+ l/2Ml/2,l = -j / / / («> X l) 

J v z;k J v yJ Jv x -i 

fi,j,k,i+i/2 = f (v x;i , v y .j, v z . k , x) dx, 

J Xl 

_ m I f X l + l r v z:k + l rVy;j+l rv x . i+ i 

M i+1/2j . +1/2)fc+1/2jJ+1/2 = — / / / / v m f {v,x)dvdx, 



-Jv z . k 



where the subscripts i,j,k, and I denote the grid position in the v x ,v y ,v z , 
and x directions, M m = (M™, M™, M™), M° x = M° y = M° = M°, M™ = 
(M™,M™,M™), Ml = Ml = Ml = M°, and m = 0,1,2. We split the 
Vlasov equation (I3~T|) into two equations in three-dimensional velocity and 
one-dimensional configuration spaces, which are alterna t ely ad vanced by the 



MMA scheme and the CIP-CSL2 scheme (jYabe et al.l . 12001), respectively . 



The Maxwell equation (J32l) is solved by the implicit scheme (jHoshinol . 11986 



as 



19871). The time integratio n of the system is carried out in the same manner 



Minoshima et al.l (120111 ). Physical variables in the system are normalized 
as follows; velocity by the speed of light, time by the inverse electron plasma 
frequency u pe , electromagnetic fields by an ambient magnetic field strength, 
and position by the Debye length A^>- The boundary conditions are periodic 
in the configuration space and open in the velocity space where constant 
incoming fluxes are assumed while outgoing fluxes are perfectly lost. The 
simulations are executed on a generic workstation with dual Intel Xeon Quad- 
Core processors. 

4-1. Perpendicular wave propagation 

We first test the linear wave propagation perpendicular to the m agnetic 
field line, which has been previously tested in Minoshima et al. ( 201lh (in the 



paper, we assumed two dimensionality in velocity space). Since the three- 
dimensional scheme is not designed in the same way as the two-dimensional 
one, we test the same problem again. The initial plasma condition is a 
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uniform and isotropic Maxwell distribution with a small (1%) uniform ran- 
dom perturbation only for the electron density. A uniform magnetic field 
is initially imposed in the z-direction. The initial electric field is deter- 
mined from the Poisson equation (Gauss's law). Simulation parameters 
are as follows; a mass ratio m p /m e = 16, a ratio of the electron gyro to 
plasma frequency uj ge /u pe = 0.5, and electron and proton thermal velocities 
Ve-,th = 0.1,f p; ^ = 0.025, corresponding to electron and proton plasma beta 
values being (3 e = /3 P = 0.04. The simulation domain in the velocity space 
is [— Avth, 4utfc] with 32 grid points in each direction for each species. The 
grid size in the configuration space is equal to Xd, and the spatial length is 
512A D . The time step is 0.05/V2O;- 1 . 

Fig. M shows the Fourier spectrum of the electrostatic field E x integrated 
until (a) uj pe t = 361.3 and (b) u pe t = 723.4. Similar to the previous sim- 
ulation, we can clearly identify the electron and ion cyclotron (Bernstein) 
modes, X- and Z-modes, and lower-hybrid waves. During the simulation 
(electrons gyrate more than fifty times), the total energy is conserved within 
an error of 0.1%. 

4-2. Parallel wave propagation 

We next test the linear wave propagation parallel to the magnetic field 
line. The initial plasma condition is a uniform and isotropic Maxwell distri- 
bution. A uniform magnetic field is initially imposed in the ^-direction, and 
then a small (1%) uniform random perturbation is added to the transverse 
field. The initial electric field is zero. Simulation parameters are the same as 
in Section [4.11 The simulation domain in the velocity space is [—4^,4^] 
with 32 grid points in each direction for each species. The grid size in the 
configuration space is 4A.D, and the spatial length is 2048A£>. The time step 

is O.l/VHe 1 . 

Fig. 17(a) shows the Fourier spectrum of the transverse field B y integrated 
until u pe t = 723.4. For comparison, we also perform the electromagnetic PIC 
simulation with the same parameters (except that the grid size is Xd in the 
PIC), and the result is shown in Fig. [TJb). The number of particles in each 
cell is 12,500 so that the total memory usage is comparable between the 
two simulations. We can clearly identify the R- and L-modes, and whistler 
waves. The ion-cyclotron wave is not clear because the integration time is not 
sufficiently long. During the Vlasov simulation, the total energy is conserved 
within an error of 0.005%. 
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The high frequency whistler waves (cu > —kv e . >t h+u ge ) effectively dissipate 
in the Vlasov simulation through the cyclotron damping by thermal electrons, 
while it is not clear in the PIC simulation owing to the thermal noise. 

4-3. Electron temperature anisotropy instability 

We lastly test the nonlinear evolution of whistler w a ves t hrough the 



electron temperature anisotropy instability (jSydora et al.l . 120071 ). The ini- 
tial condition is a uniform and isotropic Maxwell distribution for protons, 
and bi-Maxwell distribution for electrons with a temperature anisotropy 
T e ±/T e \\ > 1, where T e ± and T e \\ are temperatures perpendicular and par- 
allel to the magnetic field line. A uniform magnetic field is initially imposed 
in the x-direction, and then a uniform random perturbation is added to 
the transverse field to initiate the instability. The initial electric field is 



zero. Simulation parameters are as follows (same as in lSydora et al.l ( 120071 )): 
m p /m e = 1836, LO ge /u pe = 0.2, /3 e n = f} p — 1, and T e ±/T e \\ = 3. The actual 
mass ratio is employed because protons do not play an important role in this 
instability. The simulation domain in the velocity space is [— 4.5?; e; ^, 4.5i> e;i /J 
for electrons, and [—3v p;t h, 3v P]t h\ for protons with 32 grid points in each di- 
rection. The grid size in the configuration space is Ax = 4Ad, and the spatial 
length is L = 2048A D . The time step is 0.25/^/2^. 

Fig. E^a,b) shows the time profile of the transverse electric field spectrum 
E z (k,t) and distribution E z (x,t). At the linear phase, we observe wide-band 
waves in the wavenumber range of kc/u pe = 0.5 — 1.0. During the nonlinear 
phase, the wavelength shifts to longer one (kc/u pe = 0.3 — 0.4), and nearly 
coherent waves propagate forward a nd backward. These features are in good 
agreement with Sydora et al. ( 20071 ). Fig. EJc) shows the Fourier spectrum 



of the transverse electric field E z (k, u) superimposed on the linear dispersion 
relation of the whistler wave (dashed line). The excited waves are certainly 
the whistler waves. Fig. E£d) shows the time profile of E z at the wavenumber 
corresponding to the fastest g rowing mode (kc/u pP = 0.7) . The growth rate 



agrees with the linear theory (IGurnett and Bhattacharjed (120051 )). 



Fig. [H^a,b) shows the longitudinal electron distribution function ff f e dv y dv z 
at linear and nonlinear phases. The instability increases the electron temper- 
ature parallel to the ambient magnetic field line. Fig. M^c) shows the time 
profile of the spatially-averaged perpendicular temperature T e ± = (T ey + 
T ez )/2, parallel temperature T e \\ = T ex , and temperature anisotropy T e ±/T e \\. 
The temperature anisotropy is decreased as the electric field is increased (see, 
Fig. Eld)). At the nonlinear phase, the system reaches marginal stability. 
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The saturation level of the tempera ture anisotropy T e x/T e \\ ~ 1.2 is in good 



agreement with ISydora et al.l (120071 ) . Fig. WA) shows the time profile of the 
total energy with different spatial resolution (Ax, L). A black line shows the 
reference result (4Ad, 2048A£i), a red line with finer grid size and narrower 
length (\e>, 512Xd), and a blue line with narrower length (4Ad, 512Ad). The 
total energy is slightly increased by 1.6%. The error level is not dependent 
on the spatial resolution. Therefore, we consider that the error is mainly 
caused by an inaccuracy of the distribution function in velocity space. 

5. Summary and discussion 

We have presented an extens ion of the multi-moment advection (MMA) 



scheme ( jMinoshima et al.l . 120111 ) to the three-dimensional case, for full elec- 
tromagnetic Vlasov simulations of magnetized plasma. The scheme treats 
not only point values of a profile but also its zeroth to second order piecewise 
moments as dependent variables, and advances them on the basis of their 
governing equations. Similar to the one- and two-dimensional schemes, the 
three-dimensional scheme has quite high capability for Vlasov simulations. 

The scheme is applied to the linear and nonlinear electromagnetic Vlasov 
simulations. Since the scheme can solve the solid body rotation and ad- 
vection problem with little numerical dispersion or diffusion, it enables us 
to perform long time Vlasov simulations of magnetized plasma with small 
numerical errors. In the electron temperature anisotropy instability (Sec- 
tion I4.3p . our Vlasov simulation code successfully describes the cooling as 
well as heating processes and the marginally stable state, by virtue of the 
diffusionless property of the scheme. 

In this scheme, we apply the Weighted ENO construction technique for 
line-integrated variables (eq. f fT9|) ). This does not mean that the scheme 
possesses the non-oscillatory property. Better techniques may be devised to 
suppress numerical oscillations. 

In the solid body rotation problem (Section 12. 2p . we split the velocity into 
two phases. Since the order of the first and second phases is arbitrary, a sim- 
ulation result is not necessarily same when one reverses the order. However, 
we confirm tha t the effect is negl i gible small. 



As shown in lMinoshima et al.l ( 120111 ). the one- and two-dimensional MMA 



schemes exactly guarantee the conservation of the zeroth to second order cen- 
tral moments in the advection problem with constant velocity, and the con- 
servation of the sum of the second order moments in the solid body rotation 
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problem. The same holds for the tree-dimensional scheme in the advection 
problem, however, not in the solid body rotation problem, owing to the split- 
ting procedure. Nevertheless, benchmark tests have shown that the scheme 
preserves the profile and the orbit of rotation with high accuracy. 

One of advantages of the Vlasov simulation against the (explicit) PIC 
simulation is that the grid size in configuration space is not necessarily re- 
stricted to the Debye length. In fact, we set the grid size larger than the 
Debye length in Sections 14.21 and 14.31 This advantage can be applied espe- 
cially to the simulation with the large frequency ratio (cu pe /cu ge 3> 1). Due to 
the restriction of the grid size, the frequency ratio u pe /u ge in many explicit 
PIC simulations is much smaller than in our space environment, to save com- 
putational cost. The Vlasov simulation can be performed with larger u pe /u ge 
within reasonable computational cost by using a coarser grid, unless Debye- 
scale structures are important. Therefore, our Vlasov simulation code will 
be able to simulate large-scale and long-time plasma kinetic phenomena with 
large u pe /uj ge . Another advantage of the Vlasov simulation is the simplicity 
for parallel computation, because both the plasma and electromagnetic fields 
are treated as Eulerian variables. In these points of view, the Vlasov simu- 
lation is a necessary technique for the plasma kinetic simulation on present 
peta-scale and future exa-scale supercomputer systems. 
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Appendix A. Coefficients of the interpolation function of MMA3D 

Clll;i,j,k = 
Cll2-iJ,k = 

C\li;i,j,k = 
Cl22;i,j,k = 
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Remaining coefficients can be obtained on the basis of a cyclic rule. For 
example, replacing (x,y,z) — > (z,x,y), (£, 77, 0) — >■ (8,CiV)i an d {h3>fy ~~ 
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(k,i,j) in Equation (|A.5p (e.g., l x - iC eu,j up ,k ->■ lz;iu P ,j,kceii) S ives Cai^ij,*- 
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(a) Initial (b) MMA (c) CSL2 




Figure 1: Three-dimensional solid body rotation problem of a symmetric gaussian profile, 
(a) Initial profile. (b,c) Profiles after 100 rotations calculated with the MMA and CIP- 
CSL2 schemes. 
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Figure 2: Errors of the three-dimensional solid body rotation problem of a symmetric 
gaussian profile as a function of the grid size. Triangles and diamonds are obtained from 
the MMA and CIP-CSL2 schemes. A dashed line indicates the third order accuracy. 
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(a) Initial 



(b) MMA 



(c) CSL2 



Figure 3: Three-dimensional solid body rotation problem of an asymmetric gaussian pro- 
file. The format is same as Fig. [TJ 
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Figure 4: Temporal variation of the standard deviation (a) <j x , (b) a y , (c) a z , and (d) their 
average in the three-dimensional solid body rotation problem of an asymmetric gaussian 
profile. Solid and dashed lines are obtained from the MMA and CIP-CSL2 schemes. 
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Figure 6: Electromagnetic Vlasov simulation of perpendicular-propagating waves. Color 
contour shows the Fourier amplitude of the electrostatic field E x (normalized by its max- 
imum value), (a) The simulation result until u pe t = 361.3. Horizontal and vertical axes 
are the wavenumber and frequency normalized by the inverse electron gyro radius and the 
electron gyro frequency. The amplitude is exponentiated by 0.15 for illustration. From top 
to bottom, dashed lines represent the R-mode cutoff, upper hybrid, L-mode cutoff, and 
lower hybrid frequencies. A dot-dashed line represents the dispersion relation of the light 
mode in vacuum, (b) The simulation result until oj pe t = 723.4. Horizontal and vertical 
axes are the wavenumber and frequency normalized by the inverse proton gyro radius and 
the proton gyro frequency. The amplitude is exponentiated by 0.3. A dot-dashed line 
represents the dispersion relation of the Alfven wave. 
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Figure 7: Electromagnetic (a) Vlasov and (b) PIC simulations of parallel-propagating 
waves. Color contour shows the Fourier amplitude of the transverse field B y (normalized 
by its maximum value). Horizontal and vertical axes are the wavenumber and frequency 
normalized by the inverse electron gyro radius and the electron gyro frequency. From 
top to bottom, dashed lines represent the R-mode cutoff, L-mode cutoff, and electron 
gyro frequencies. Dot-dashed lines represent the dispersion relation of the light mode in 
vacuum and the Alfven wave. 
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Figure 8: Electromagnetic Vlasov simulation of the electron temperature anisotropy in- 
stability. (a,b) Time profile of the transverse electric field spectrum E z {k,t) (normalized 
by its maximum value) and distribution E z (x,t). (c) Fourier spectrum of the transverse 
electric field E z (k,uj) (normalized by its maximum value). A dashed line represents the 
linear dispersion relation of the whistler wave, (d) Time profile of E z at the wavenumber 
corresponding to the fastest growing mode (kc/tu pe = 0.7). A dashed line indicates the 
linear growth rate. 
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Figure 9: Electromagnetic Vlasov simulation of the electron temperature anisotropy insta- 
bility. (a,b) Longitudinal electron distribution function at linear (uj ge t = 28.3) and nonlin- 
ear (uj ge t = 283) phases, (c) Time profile of the perpendicular temperature (dashed line), 
parallel temperature (dot-dashed line), and temperature anisotropy (solid line). The tem- 
peratures are normalized by the initial value of the parallel temperature, (d) Time profile 
of the total energy (normalized by its initial value). Black, red, and blue lines are obtained 
from different simulation runs with (Ax/Xd,L/X d ) = (4,2048), (1,512), and (4,512), re- 
spectively. 
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